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ABSTRACT 

We study the fate of internal gravity waves approaching the centre of an initiahy 
non-rotating solar-type star, by performing three-dimensional numerical simulations 
using a Boussinesq-type model. These waves are excited at the top of the radiation 
zone by the tidal forcing of a short-period planet on a circular, coplanar orbit. This 
extends previous work done in two dimensions by Barker & Ogilvie. We first derive 
a linear wave solution, which is not exact in 3D; however, the reflection of ingoing 
waves from the centre is close to perfect for moderate amplitude waves. Waves with 
sufficient amplitude to cause isentropic overturning break, and deposit their angular 
momentum near the centre. This forms a critical layer, at which the angular velocity 
of the flow matches the orbital angular frequency of the planet. This efficiently absorbs 
ingoing waves, and spins up the star from the inside out, while the planet spirals into 
the star. 

We also perform numerical integrations to determine the linearised adiabatic tidal 
response throughout the star, in a wide range of solar-type stellar models with masses 
in the range 0.5 ^ tti^/Mq ^ 1.1, throughout their main sequence lifetimes. The 
aim is to study the influence of the launching region for these waves at the top of 
the radiation zone in more detail, and to determine the accuracy of a semi-analytic 
approximation for the tidal torque on the star, that was derived under the assumption 
that all ingoing wave angular momentum is absorbed in a critical layer. 

The main conclusions of this work are that this nonlinear mechanism of tidal 
dissipation could provide an explanation for the survival of all short-period extrasolar 
planets observed around FGK stars, while it predicts the destruction of more massive 
planets. This work provides further support for the model outlined in a previous 
paper by Barker & Ogilvie, and makes predictions that will be tested by ongoing 
observational studies, such as WASP and Kepler. 

Key words: planetary systems - stars: rotation - binaries: close - hydrodynamics - 
waves - instabilities 



1 INTRODUCTION 

Tidal interactions are important in determining the fate of 
short-period extrasolar planets and the spins of their host 
stars. The extent of the spin-orbit evolution that results from 
tides depends on the dissipative properties of the star and 
planet. These are usually parametrized by a dimensionless 
quality factor for each body, which is an inverse measure of 
the dissipation. This is usually defined to be proportional to 
the ratio of the maximum energy stored in a tidal oscillation 



to the energy dissipated over one cycle (e.g. ^Goldreich 
Soter|1966| ). 



In principal the modified quality factoi[^(3^ depends on 
tidal frequency, the internal structure of the body, and the 
amplitude of the tidal forcing. Unfortunately, the mecha- 
nisms of dissipation that contribute to Q' are poorly under- 
stood. To simplify this difficult problem, nearly all studies 
neglect any amplitude dependence of the dissipation (ex- 
cept, for example Goodman Lackner|[2009 ). Such studies 
already exhibit a complicated dependence of Q' on the tidal 



frequency (e.g Savonije Papaloizou 1997jWu|2005 Ogilvie 
|fc Linll2004l hereafter OL04; [Ogilvie Lin||2007p iereafter 
OL07; |Papaloizou Ivanov|2010| , and the internal structure 
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^ Related to the traditional Q by Q' = 3Q/2k, where k is the 
second-order potential Love number of the body. 
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of the body (e.g 0LQ7; [Barker Ogilvie|2009|). In this paper 
we extend the work of Barker & Ogilvie (2010|) (hereafter 
BO 10) to study an important nonUnear dissipation mecha- 
nism in solar- type stars. 

We can decompose the tidal response of a fluid body, 
which in this paper we take to mean a solar- type star, into 
an equilibrium and a dynamical tide, defined such that the 
total displacement is the sum of these two displacements, i.e, 
that ^ = C'^ + C^- The equilibrium tide is a quasi- hydrostatic 
bulge defined by 



5 

g 



and V-C 



0, 



(1) 



in stratified regions dGoldreich Nicholson] 1989| ), where ^ 
is the tidal gravitational potential experienced by the body 
and g is the gravitational acceleration. The total displace- 
ment is not well described by Eqs. [l] in convective regions 
(Goodman & Dickson 1998, hereafter GD98; |Terquem et al 



1998[ hereafter T98; 0L04). This is because in a barotropic 
flow (with adiabatic stratification) vorticity is conserved, so 
we must have V x ^ = 0, whereas V x 7^ 0, in general. 
The presence of a convection zone (hereafter CZ) thus im- 
plies that a dynamical tide must exist. The dynamical tide 
is defined as the residual response that results from the 
equilibrium tide not being the exact (linearised) solution to 
the problem, when the tidal frequency is nonzero. 

The equations governing the adiabatic equilibrium and 
dynamical tides in linear theory are 



= 



-UJ I 



(2) 
(3) 



from which it is clear that the dynamical tide is not forced 
directly by the tidal potential, only by the inertial terms in 
the equation of motion. 

Dynamical tides in radiation zones of solar-type stars 
take the form of internal (inertia-) gravity waves (IGWs), 
which have frequencies below the buoyancy (or Brunt- 
Vaisala) frequency N. These have previously b een proposed 
for early- type stars (e.g. Zahn 1975), 



to contribute to 
where they are damped at the surface by radiative diffu- 
sion. It has also been proposed that these waves could syn- 
chronise the spin of the star with the orbit (in this case of 
a close-binary perturber) from the outside in (Goldreich & 
Nicholson|1989t . In BOlO, we considered a nonlinear mech- 
anism of tidal dissipation in solar-type stars (with radiative 
cores), extending an idea by GD98. The consequences of this 
mechanism are similar to [Goldreich & Nicholson (1989), ex- 
cept that the star is synchronised with the orbit from the 
inside out. 

A short-period planet excites IGWs at the top of the 
RZ of a solar- type star, where N increases from zero with 
distance into the RZ from the CZ/RZ interface. There is 
thus a location at which N ^ 1/^, with P being the plane- 
tary orbital period, at which IGWs (which have frequencies 
less than N) are efficiently excited. These waves propagate 
downwards into the radiation zone (hereafter RZ), until they 
reach the centre of the star, where they are geometrically fo- 
cused and can become nonlinear. If their amplitudes are suf- 
ficiently large, convective overturning occurs, and the wave 



breaks. This has consequences for the tidal torque, and the 
stellar Q^ In this paper we study this mechanism, primarily 
using three-dimensional numerical simulations. 

In BO 10, we derived a Boussinesq-type system of equa- 
tions that is relevant for describing the dynamics of IGWs 
approaching the centre of a solar- type star. We then per- 
formed numerical simulations, solving these equations in two 
dimensional cylindrical geometry. In this paper we extend 
these simulations to three dimensions, in spherical geometry, 
and confirm that the most important results of BO 10 are not 
affected by this extension. We first derive a linear wave solu- 
tion, which represents the waves excited by planets on short- 
period orbits as they approach the central < 5% of the star, 
within which N oc r and g (x r. A weakly nonlinear analysis 
confirms that this solution is not exact, though the reflection 
of ingoing waves is close to perfect for moderate amplitude 
waves. We present the numerical setup and the analysis of 
the simulation results. We then discuss the launching region 
at the top of the RZ, and the possible effects of magnetic 
fields on the problem. Together with BO 10, we provide a 
possible explanation for the survival of all short-period ex- 
trasolar planets observed thus far, which will be tested by 
ongoing observational studies, such as WASP and Kepler. 



2 TIDAL POTENTIAL 

The tidal potential experienced by a star hosting a short- 
period planet can be written in standard spherical polar 
coordinates (r, 0,(j)) as a sum of rigidly rotating spherical 
harmonics 



Re 



r Yi 



(4) 



in a non-rotating (but non-inertial) reference frame centred 
on the star, where is a spherical harmonic (normalised 
such that the integral of lY^^I^ over solid angles is unity) and 
^i,rn is an amplitude. Here cj is the frequency in that frame, 
related to the tidal frequency by a; = cj — mQ, where Q is the 
spin angular frequency of the star. In this paper we consider 
the waves excited by planets on circular, coplanar orbits, 
relegating any studies of the waves excited by eccentric and 
inclined planets to future work. In this case the dominant 
component of the tidal potential is quadrupolar (/ = 2), with 
m — 2, and takes the form 



^(r,6',(/),t) 



(5) 



where rup is the planetary mass and a is its orbital semi- 
major axis. Since most short-period planets orbit faster than 
their stars spin, as a result of stellar magnetic braking, we 
take the star to be (initially) non-rotating, so that cj = cj = 
2n, where n is the orbital angular frequency of the planet. 
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3 LINEAR THEORY 

In this section we derive a linear wave solution, starting from 
the adiabatic Boussinesq-type system (BO 10) 



Du= -\/q + rb, 
Db + C'^r -u = 0, 
V-u = 0, 



(6) 

(7) 

(8) 
(9) 



where u is the fluid velocity, 6 is a buoyancy variable (pro- 
portional to the entropy perturbation) and ^ is a modified 
pressure variable. This system of equation was derived in 
2D but is equally valid in 3D. We note that in this model 
N — Cr^ where C is a constant that measures the strength 
of the stable stratification at the centre. This model is valid 
in the innermost < 3% of the star, which contains multiple 
wavelengths for the gravity waves excited by short-period 
planets. We define — uj/m^ which is the pattern speed of 
the forcing, and equals the orbital angular frequency of the 
planet. We furthermore use the non-dimensionalisation that 
the unit of length is Vi^jC and the unit of time is XjVL^. 
Linearising about hydrostatic equilibrium in 3D spherical 
geometry, we seek solutions steady in the frame rotating at 
the angular rate Qp, proportional to e*^^, where ^ = (\) — VL^t 
is the azimuthal coordinate in this frame. This leads to the 
following equations: 



-imur 



-drq + rb, 
1 , 
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-imb — 



r smt/ 



1 2 1 

^dr(r Ur) H ^ 

rsmt 



(10) 

(11) 

(12) 
(13) 

: 0. (14) 



To obtain the linear solution, we expand scalar quantities 
(i.e. q and b) in spherical harmonics Y{^{0^ ^), since the prob- 
lem is separable in both angular coordinates. In Eqs. p^}fT4l 
we have already included the ^-dependence of these func- 
tions (e'"^^). We thus take q = q{r)Yi^{e, 0), so that the to- 
tal functions are expanded onto Yi^{0^^)^ and similarly for 
b and Ur^ where from now on we drop the hats on the radial 
functions. The remaining velocity components are expanded 
onto angular functions as appropriate to satisfy Eqs. [To}{l4] 
The relation 



-dr(r'^Ur), 



" /(/ + !) 
follows from incompressibility and the result 

1 . . ..m. . 1 



(15) 



sm tf sm 
This enables us to derive the linear differential equation 



(16) 



r )Ur 



0, 



(17) 



whose solutions can be written in terms of Bessel functions 
of half-integer order (alternatively these can be written as 
spherical Bessel functions, or they can be reduced to elemen- 
tary functions) . The corresponding (total) linear solution for 
a standing wave in 3D can be written (where real parts are 



assumed to be taken) 
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(18) 



^-9. [r4j,+ i (fcr)]a,yr(^,0, (19) 
{kr)Yr{e,i), (21) 



where B G C is an amplitude, and 
-V/(/ + l). 



(22) 



Ingoing wave (hereafter IW) and outgoing wave (hereafter 
OW) solutions can be obtained by replacing the Bessel func- 
tion of the first kind by equivalent Hankel functions of the 
first [Ju + lYu) and second kinds [Ju — iYu), respectively. 
Note also that 



r 



r 2 



[(1 + i (fcr)±iy,+ i (fcr)) 
kr (kr) ± (fer))] . (23) 



Starting from Eq. 46 of BO 10, we can calculate a con- 
served energy flux. Integrating this equation over ^ elim- 
inates the terms containing derivatives in t and due to 
periodicity in ^, as a result of the fundamental theorem of 
calculus. This allows the definition of a conserved quantity 
proportional to the energy flux. 



F= / sin Our{E + q)d^dO, 
Jo Jo 



(24) 



where the energy density per unit volume is ^ = 
|po + with po the (constant) central density of 

the star. For linear waves, terms involving products are 
small, so ^ ^ ^, to this order. This leaves 



F — Tvr' 



/ Re [urq*] 
Jo 



sin OdO. 



(25) 



Whether this is positive or negative depends on whether the 
wave is ingoing or outgoing. 

Substituting the linear solution into Eq. [25] provides a 
simple expression for the flux of a single /, m wave: 



F= . A \Ao^tf -\Ai„\'), 



TVIH + I) 



(26) 



with corresponding energy flux — p^F ^ and angular mo- 
mentum flux F^ — P^F. We define the complex ampli- 
tudes of the IW and OW to be Ain and Aout^ respectively. 
Eq.|26| follows from the orthonormality of spherical harmon- 
ics and the Wronskian of the Hankel functions 



Jo Jo 

W[J^.{kr) + iY^^ikr), J^{kr) - iY^{kr)] = 



sinOdOd^^dW^ . (27) 



(28) 



The particular wave that we will study has / = m = 2, and 
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the components 
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--a. [rij5(fcr)]aey2'(^,0, (30) 



(32) 



where k — 



3.1 Criterion for overturning isentropes 

A condition in 3D for isentropic overturning can be derived 
from considering when the radial gradient of the entropy 
s = 6+ (l/2)r^ becomes negative, i.e., when (l/r)drh < —1. 
This is equivalent to the criterion that \{l/r)dr{rur)\ > \m\, 
which can be shown by using Eq. |13| To correlate our nota- 
tion with the appendix of OL07, we define the dimensionless 
nonlinearity parameter A, such that overturning occurs if 
\A\ > 1. This is defined such that the radial velocity is the 
real part of 



40Ar~ 
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7E 



1 



-r^ 1 sin/cr- 
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-r cos kr 



sin 9e 



(33) 



in the dimensionless units that we have been using in this 
section (this is equivalent to Eq. 



29) 



Overturning is achieved at the centre when the radial 
velocity in the wave takes a maximum value Ur — 1.27 at its 
innermost peak at r = 2.04. This value is used to compare it 
to the magnitude of the radial velocity achieved in numerical 
simulations at the onset of wave breaking. In these dimen- 
sionless units, we similarly require > 0.99 or b > 0.38. 
Note that there is not such a simple interpretation of the cri- 
terion on U(f) as in 2D, though the value is quantitatively very 
similaij^ From the results of the 2D simulations in BO 10, 
we expect the waves to undergo instability and break within 
several wave periods after these criteria begin to be satisfied. 



3.2 Weakly nonlinear theorjj^ 

The linear solution written down in Eqs. |18f|21| is not a 
nonlinear solution, unlike the equivalent in 2D. This can 
be shown by computing the nonlinear terms in the full 
Boussinesq-type system using the linear solution, in Mathe- 
matica, for example. We find that u- V6 ^ 0, in general, and 
similarly for the nonlinear terms in the momentum equation. 
This means that the reflection of the waves from the centre 
could be different than in 2D, since nonlinearities do not 
vanish for this wave. In this section we perform a weakly 



^ In the figures below, we use a different nondimensionalisation of 
the velocity, by normalising it to the constant (asymptotic) radial 
phase velocity uoXr/ 2tt in the wave, this is identical to BOlO, and 
gives values exactly 1/2 of those in the units of this section. 
^ I would like to thank Gordon Ogilvie for suggesting the calcu- 
lations in this section. 



nonlinear analysis to determine the dominant nonlinear ef- 
fects for small amplitudes. Since these nonlinearities do not 
vanish, this highlights the importance of numerical simula- 
tions for these waves approaching the centre. We describe 
the results of such simulations in § [5] 

We propose a weakly nonlinear solution of the form 

+e^{ur2o{r, 6) 

^ (wr22(r, 6')e^*^ + ul22{r, 6')e"^'^) } 
+0(e^), 



+ 



(34) 



and similarly for the other variables, where e ^ 1. This 
form is adopted because we are interested in calculating 
whether the incoming wave generates harmonics through 
the quadratic (self-) nonlinearities. These additional waves 
(other than m = 0) will escape to infinity and carry away 
a portion of the energy flux. Here we write Uri{r^O)e^^ — 
Ur{r, 0,^) from Eq. |29] for the I — m — 2 wave above, and 
similarly for other variables. We substitute these expansions 
into the Boussinesq-type system and equate powers of e. At 
each order we also equate coefficients of e*^^ . At leading or- 
der only one mode is present, and we obtain the previously 
derived linear solution. After some algebra the solution at 
O(e^) can be computed to give 



A 



^ [J9/2 (y^kr) 
+il9/2(y576^r)]y4'(^,0, 



(35) 



which is an / = m = 4 wave with complex amplitude A22, 
which has been computed using Mathematica and given in 
terms of A. 

For the wave described by Eq. [35] F can be computed. 
The ratio of the energy flux in the outgoing I = m — A wave 
to the ingoing I — m — 2 wave can be shown to be approxi- 
mately 1.2 X 10~^|A|^. We deflne a reflection coefficient 



7^: 



Ai'^ 



(36) 



which measures the amplitude decay for a wave travelling 
from a radius r to the centre, and back to r. For perfect re- 
flection from the centre IZ — 1^ whereas complete absorption 
means that IZ — ^. The reflection coefficient for reflection 
from the centre, for a weakly nonlinear I — m — 2 wave, can 
be computed from 



1.2 X 10"^|Ap. 



(37) 



This means that a weakly nonlinear primary wave (with 
\A\ <^ 1), will reflect approximately perfectly from the cen- 
tre, with a reflection coefficient that is close to unity. How- 
ever, a small fraction of the IW energy flux is transferred 
to waves with higher / and m- values, reinforcing the fact 
that Eqs. |18f|21| is not an exact solution, contrary to the 
analogous solution in 2D. 



4 NUMERICAL SETUP 

We solve the Boussinesq-type system in three dimensions us- 



ing the Cartesian pseudospectral code SNOOPY (Lesur & 



Wave breaking in a solar-type star 5 



Longaretti|2005| ), as in BO 10 (see § 6 of that paper for fur- 
ther details). However, we modify the forcing and damping 
to take into account the z-direction, and instead of forcing 
an m = 2 wave in the equatorial plane, we now have 



f = -frRe[Y2^{e,^-'^t)]er, 



y^) COS (jjt — 2xysm(jjt} e^, (38) 

in Cartesian coordinates (with — + + z^), which 
is applied in the region O.'^brhox ^ t ^ O.Orbox- We study a 
region x,y,z G [-rhox.rhox], where rhox = 1-5, in arbitrary 
units (not the same as §[3]). For r > O.^Vhox the solution is 
damped to zero by using a parabolic smoothing function, as 
in BOIO. 

We primarily use a resolution of 256^, for which the 
simulations were possible to run on a single Intel Core 17 
machine, utilising all 8 cores, with a typical run time of sev- 
eral weeks to resolve a hundred wave crossing times. We set 
uj — 1 and choose a typical IGW wavelength of Ar = 0.15, 
giving approximately 8 wavelengths within the box. The 
value of Ar is chosen to be slightly larger than that used 
in most of the 2D calculations. This increases the number 
of grid points within a wavelenth of the primary wave, par- 
tially offsetting the reduction in resolution that results from 
using a smaller number of grid points per dimension than 
in 2D. Choosing larger wavelengths than this is found to 
result in unwanted effects from the proximity of the forcing 
region, which modifies the linear solution. The viscosity and 
radiative diffusion coefficients, at least one of which must be 
implemented in the code for numerical stability, are chosen 
to he V — Ax 10~^ and = 0. An otherwise identical setup 
is used as in BOIO, except using spherical geometry instead 
of cylindrical geometry. We have confirmed that the linear 
solution is well reproduced with this numerical setup. 



5 NUMERICAL RESULTS 

In this section we describe the results of the numerical 
simulations. We analyse the results of the simulations us- 
ing the IW/OW decomposition described in Appendix [A] 
and reconstruct the solutions from the computed IW/OW 
amplitudes to compare with the simulation data, as de- 
scribed therein. From preliminary investigation, we find that 
fr = fr27v/{(jo'^Xr) > 0.2 Is rcqulrcd for breaking, so a variety 
of simulations are performed with forcing amplitudes either 
side of this value. The basic results of these simulations are 
that the wave reflects approximately perfectly from the cen- 
tre of the star if the amplitude of the wave is smaller than a 
certain critical value, which is found to correspond with that 
required for isentropic overturning. Above this value, wave 
breaking and critical layer formation occur. This picture is 
identical to that in 2D. 



5.1 Low-amplitude simulations 

For a low-amplitude simulation, with max(itr) below the 
critical value for isentropic overturning (using fr = 0.1), we 
plot the variation in amplitudes of the IWs and OWs, and 
also the reconstructed solutions in Fig ^ We analyse the 
results using the method described in Appendix [A] choosing 
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Figure 2. 2D plots of Ur (top) and (bottom) on the xy-plane 
for a simulation in which the wave perfectly reflects from the 
centre, for which fr = 0.1 This can be qualitatively compared 
with Fig. 3 in BOIO. 



a time once transients have been sufficiently damped and 
standing waves have formed. The decay in radius is roughly 
(though slightly smaller than) that which would be expected 
from viscous damping, by the fractional amount 



Ur 
Ur,0 



: exp 



Jo ^9,r 



dr 



(39) 



where c^,r = CA^/(27r^) is the (constant asymptotic) radial 
group velocity, for a wave of the given wavelength. As in 
2D, we have confirmed this explanation by running simula- 
tions without viscosity, which are found to not exhibit this 
decay (though these simulations eventually become numeri- 
cally unstable liv — n — Qi). Using a smaller viscosity is also 
found to reduce the wavelength-scale oscillations around the 
mean slope. These result from the fact that the linear so- 
lution to the forced wave problem is no longer exact in the 
presence of viscosity. We find that increasing the number 
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Figure 1. In the top left panel we plot Ain (blue) and Aout (red) vs. r, from a low- amplitude simulation in which the primary wave 
approximately perfectly reflects from the centre, with IZ displayed in the top right panel. Below, are the velocity components and the 
buoyancy variable (blue), plotted together with the reconstructed linear solution (red), against radius. This is taken from a simulation 
with fr = 0.1 at t = 250, once standing waves have formed. The amplitude decay with radius can be explained as due to viscous damping. 



of grid points within each shell (by reducing the values of 
i step step and kstep) slightly reduces the vertical extent of 
these oscillations, because this averages out the errors that 
result from the assumption that the inviscid linear solution is 
exact. However, increasing the number of grid points within 
each shell has negligible effect on the mean slope. 

The spatial structure of the solutions in three dimen- 
sions in the x?/-plane is very similar to that in two dimen- 
sions, as can be seen in Fig. [2] (which can be compared with 
Fig. 3 in BOlO). In Fig.|3] we plot tt^ on the x^-plane. This 
shows that the magnitude of the azimuthal velocity peaks 
at — 7r/2, due to the latitudinal form of . 

From the calculation in § |3.2| we expect the effects of 
nonlinearity to be much weaker than the effects of viscos- 
ity for small-amplitude waves which do not cause convective 
overturning. Since the effects of weak nonlinearity are very 
small, it is difficult to quantitatively confirm the results in 
§ |3.2| using, for example, an extension of the method de- 
scribed in Appendix^ for multiple / and m values. Never- 
theless, we have qualitatively confirmed the result that the 
reflection is coherent and nearly perfect (in that 7^ ~ 1) for 
amplitudes below that required for overturning the strati- 
fication. As in 2D we do not observe any instabilities that 
act on the waves when they have insufficient amplitude to 
overturn the stratification. In this case, the waves can form 
global modes in the RZ. 



5.2 High-amplitude simulations 

In high- amplitude simulations, in which the wave amplitude 
exceeds the overturning criterion, the wave overturns the 
stratification during part of its cycle and a rapid instability 
acts on the wave, which leads to wave breaking within 1 — 3 
wave periods. This causes the rapid (within several wave 
periods) deposition of primary wave angular momentum, 
which spins up the mean flow to (which corresponds with 
the orbital angular frequency of the planet) and produces a 
critical layer. This critical layer acts as an absorbing barrier 
for IWs, as is shown from Fig. ^ which plots the variation 
in amplitude of the IW and OW, and also the reconstructed 
wave solutions (which can be contrasted with Fig. [T]). Once 
the critical layer has formed, we find |Aotit| ^ |^in|- The 
central regions are not well described by the linear model, 
as we would expect. However, the region outside oi r ^ 0.25 
is well described by the linear solution, with |Aout| ^ l^inl- 
In this region, 7^ ^ 1, so it is reasonable to assume that the 
IWs are efficiently absorbed near the centre. This picture is 
identical to that in 2D. 

The picture in 3D in the x2/-plane is very similar to that 
in 2D, as can be seen in Fig.|5](to compare with Fig. 6 from 
BOlO). However, one noticeable difference is that the pri- 
mary wave preferentially transfers its angular momentum at 
low latitudes, close to the equatorial plane. This can be seen 
in Fig. [6] where we plot the angular frequency of the flow 
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Figure 3. 2D plot of on the xz plane for a simulation in which 
the wave perfectly reflects from the centre, for which /r = 0.1. 
is of largest magnitude in the equatorial plane, where 1^2^! peaks. 



normalised to once a critical layer has formed, in both 
the xy and xz planes. This is a consequence of the latitu- 
dinal form of 1^2^ 5 whose magnitude peaks at ^ = 7r/2, as is 
illustrated in Fig.js] The critical layer absorption is observed 
to continue as the wave forcing is ongoing, so this differen- 
tial rotation is continually reinforced by the absorption of 
/ = m = 2 IWs. Since there are wave motions in the region 
of fluid interior of the critical layer, parts of these regions 
spin slightly faster than (this is not seen in Fig. |6] due to 
the adopted colour scale) . This was also observed in the 2D 
simulations. 
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6 WAVE LAUNCHING REGION AT THE TOP 
OF THE RADIATION ZONE 

The main motivation for our work is to study Q' for solar- 
type stars, and in particular to connect this with the sur- 
vival of close-in extrasolar planets. We have demonstrated 
that the fate of IGWs approaching the centre of a solar- 
type star, outlined in BO 10, is unaffected by the extension 
to three dimensions. A critical layer is formed by the de- 
position of angular momentum through wave breaking once 
the waves cause isentropic overturning near the centre. For 
lower amplitudes, the waves reflect coherently and approxi- 
mately perfectly from the centre of the star, and may form 
global modes. In that case, the dissipation is only efficient 
when the system enters a resonance with a global mode. 
The corresponding time-averaged dissipation rate is weak, 
because the system spends most of its time out of resonance 
(T98; GD98), resulting in negligible tidal evolution of the 
planetary companion. Note, however, that this neglects the 
possibility that passage through resonance will cause wave 
breaking, which should be studied in future work. 

If a critical layer forms, wave absorption is efficient, and 
global modes (of any frequency very similar to the orbital 
frequency) are prevented from being set up in the RZ. The 
tidal torque can then be computed from assuming that the 



Figure 5. 2D plot of Ur (top) and (bottom) on xy-plane for 
a simulation in which breaking occurs with = 1, at = 450. 



IWs are entirely absorbed. In this case a calculation along 
the lines of GD98 for the ingoing energy flux of the waves 
excited at the top of the RZ is required. This estimates the 
tidal torque, and thus the orbital evolution of the planetary 
companion. In this section, we perform numerical integra- 
tions of the linearised tidal response in an extensive set of 
stellar models of solar-type stars with masses in the range 
0.5 ^ m^/M© ^ 1.1, throughout their main sequence life- 
times. We aim to determine the tidal torque numerically, 
and compare it with a simple model of the launching re- 
gion at the top of the RZ, which was derived in GD98 and 
discussed in BOIO. 



6.1 Numerical computation of the linearised tidal 
response throughout the star 

In this section we solve the linearised equations governing 
the adiabatic tidal response (Eqs. [2] and [3| throughout the 
star, computing the excitation of both the equilibrium and 
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Figure 4. In the top left panel we plot Ain and Aout vs. r in a high- amplitude simulation in which the primary wave breaks (fr = 1) 
dit t = 450, with IZ in the top right panel. Below these are the velocity components and the buoyancy variable (blue), plotted together 
with the reconstructed linear solution (red), against radius. The central regions are not well described by the linear model, as expected. 



dynamical tides numerically. This allows us to determine 
the ingoing energy and angular momentum fluxes in IGWs 
launched at the top of the RZ, and to check the validity 
of approximate semi- analytic formulae for these quantities, 
presented in the next section. This is important because the 
orbital evolution of a planetary companion is determined by 
the ingoing angular momentum flux absorbed at the critical 
layer. 

We solve the following coupled ODEs for the radial and 
horizontal displacements: 



dr 



d§h 
dr 



2 dlnp 
r g dr 



/(/ + !) 



ijp'rp 



+ 



N'^ dlnp 
ruj'^ dlnp 



Tip ' 
9 



Tip 



(40) 
(41) 



An outline of the derivation of these equations is presented 
in T98. Note that we are ignoring the self-gravity of the 
entire tidal response, which is reasonable because most of 
the mass of the star is concentrated near the centre. This 
assumption is certainly valid for the dynamical tide, and is 
approximately valid for the equilibrium tide. In these equa- 
tions, we take the tidal potential in the frame rotating with 



to be equal to Eq. [H] with = <J and a; = 0, so thal|^ 



/ 



(42) 



This is the amplitude of the largest tide for a circular orbit. 
In this frame, the displacement field is separated into radial 
and horizontal (non-radial) components 



(43) 



The tidal response is further decomposed into an equilibrium 
and a dynamical tide, as defined in the introduction. 

We impose a free upper boundary, i.e., take the La- 
grangian pressure perturbation Ap = at r = i?^, so that 
the Eulerian pressure perturbation Sp = — 
relation 



^. Since the 



= ^ ( — + fr 



(44) 



follows from the non-radial equation of motion, this relates 
our variables at the surface of the star. We take an IW BC 
at the inner boundary at r ~ 0.02i^^, where we match 



^ Note that we are defining spherical harmonics in a standard 
manner, normalised so that the integral of over solid angles 

is unity, unlike T98. 
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Figure 7. Typical values of the real parts of the radial (ou^r, in 
blue solid lines) and horizontal {u:^h^ in red dashed lines) velocity 
components, in ms"-*^. We use Model S of the current Sun, and 
consider the tidal perturber to be a P = 1 d, rup = IMj, planet 
orbiting the current Sun. 
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Figure 6. 2D plot of the angular velocity U(j)/R of the central 
regions in the xy (left) and xz (right) planes, normalised to the 
angular pattern speed (i.e., orbital angular frequency), in a large- 
amplitude simulation with = 1, at t = 450. Latitudinal differ- 
ential rotation is produced by absorption of ingoing I = m = 2 
waves. 

and onto the analytic solution for an IW derived in ^ 

^,(r) = A^r-i{J5{kr)+iY5{kr)), (45) 

aw = ^r-i [3{Ji{kr)+iY5{kr)) 

-kr{Jz (kr) + iYz (/cr))] , (46) 

where k — We take to be a free parameter, so that 
these equations relate the ratio of ^ and ^h- This solution is 
quite accurate when r/R^ < 5% since is negligible in this 
region. This BC is meant to represent the IW absorption at 
a critical layer. 

We solve Eqs. [40|j4T] using data interpolated from a stel- 
lar model at points required by a 4th/5th order adaptive step 
Runge-Kutta integrator, using a cubic spline interpolation. 
In particular, the coefficients of a, ^h, in Eqs. [40|j4T] are sin- 



gular at the origin, so we first multiply these quantities by 
r before interpolating their values to the locations required 
by the ODE integrator, using the stellar model parameters. 
We then divide by r, after the interpolation. This is done to 
get the correct behaviour for small r. 

Our method of solution is a shooting method to an in- 
termediate fitting point (Pr ess et al.||1992) , which we take 
to be the CZ/RZ interface, where we enforce continuity of 
the solution. The freely specifiable initial conditions for each 
ODE integration are chosen to be ^ at the surface, and 
at the inner boundary. We use ^ = as our starting "freely 
specifiable" estimate at the surface, which is an accurate ap- 
proximation since 

The Eulerian pressure perturbation for the dynamical 
tide is 

Sp" = pLo^r [a + ^i) , (47) 

from the horizontal component of Eq. [S] The radial energy 
flux at each radius is 

F = '^Im [{Sp'r^^] , (48) 

which follows from manipulating Eqs. [2] and [3] to derive an 
energy equation. 

From Eq. [S] we derive the relation 

V • Im {5/(^'')*} = pco'lm{{^y ■ t} , (49) 

which can be averaged over solid angle to give 

Im {dr (r^Sp^i^ir) } = pr^L0^lm{{^'}rC 

+i{i + mireH}- (50) 

This is telling us that the dynamical tide is forced by the 
equilibrium tide. As part of the validation of our numerical 
code, we have confirmed that this is accurately satisfied from 
the numerical solutions computed with an ingoing BC. This 
should adequately convince ourselves that the code is able 
to accurately compute the energy flux in the dynamical tide. 
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Figure 8. Numerically computed F, normalised to the prediction 
from the analytic theory in E g. |51| (solid blue line) and GD98 
Eq. 13 (red dashed line). Eq. |51| overestimates F by less than 
10%, and is thus a good approximation of the ingoing energy 
flux. 



Figure 9. A^^ normalised to ou^y^ (solid blue line), together with 
our local approximation A^^ oc (r^ — r)-*^/^ (red dashed line). This 
approximation is reasonable over a region of size ~ 0.15Rq,. We 
also plot an arbitrarily scaled proflle of F, to compare its radial 
proflle with the proflle of A^^, for P = 1 d. 



6.2 Results 

As an illustration, we present the results of our integrations 
for a fiducial case with a P = Id, rup — IMj planet orbiting 
the current Sun (for which we use Model S, described in 
Christensen-Dalsgaard et al.|1996 ) in Fig.|7| We plot the real 
parts of Lu^r (solid blue) and u^h (dashed red) throughout 
the star, which represent typical values of the radial and 
horizontal velocity components. The radial wavelength of 
the waves decreases as the waves propagate deeper into the 
RZ, where N'^ increases. The large increase in the velocity 
amplitude near the centre is evident, as predicted from the 
linear solutions in § |3] This can be compared with a similar 
calculation in T98, displayed in their Fig. 1, which our code 
correctly reproduces for the given planetary orbital period 
when a regularity condition in applied at the centre. The 
only difference between our calculations and theirs is that 
we use an IW BC, whereas they allow the waves to perfectly 
reflect from the centre. 

In Fig. [8] we plot the ingoing energy flux, normalised 
to both the semi-analytic prediction of GD98 (red dashed 
lines), and a revised expression derived in Appendix [b] (blue 
solid lines), which will be discussed in the next section. This 
illustrates that the ingoing energy flux oscillates about its 
final asymptotic value, which it eventually approaches deep 
in the RZ. 



6.3 Semi-analytical calculation of the ingoing 
energy flux 

In this section we compare the numerically computed F with 
the semi- analytic estimate of GD98, and present a slight 
refinement which is found to be appropriate for solar-type 
stars. 

In the launching region at the top of the RZ, there is 
a location at which N"^ = at r = O.TIP©, near to 

which it follows from the dispersion relation (Eq. 4 in BOlO) 
that the radial wavenumber of gravity waves vanishes. Near 
this turning point we can approximate the solution in this 
region by assuming a functional form for . GD98 take 
N"^ oc Tb — r = X, in which case the problem in the launching 
region reduces to the solution of Airy's differential equation 



for In this case GD98 derive their Eq. 13 for the resulting 
ingoing energy flux (Eq. 92 in BO 10, which we denote herein 
as Fx). This should match the asymptotic numerical com- 
putation of F. We find that a slightly better approximation 
is to take N"^ oc ^/x over the launching regiorPj since this is 
valid over ^ O.ISP© from the interface, as we illustrate in 
Fig. |9] The slope of the curve ^/x can be obtained through 
fitting to the profile of N"^ in the stellar model. This allows 
a slightly more accurate calculation of F based on the stel- 
lar model than is obtained through direct application of Fx 
(where the gradient dN^ /dx is not uniquely defined). Nev- 
ertheless, the differences are only a factor of two at most. 
In Appendix [B] we calculate the solution in the launching 
region and the resulting P, for any power law profile of N'^ 
with a positive exponent, using the framework of OL04. Us- 
ing our approximation, the radial extent of the launching 
region Liaunch ^ O.OSRq. 

fx, we obtain, for 



For the case N'^ oc 
background. 



non-rotatmg 



F = F, 



2[r(|)] 



.[/(/ + i) 



Pbrt 



dN' 



dy/x 



dx 



(51) 



We express dx^r = ctc^/cj^^^, and compute ctc by solving 
GD98 Eq. 3 in the CZ, using a linear shooting method, with 
the boundary conditions that ^r(^6) = ^r(^0) = 0- This 
matches the numerically computed asymptotic value of P for 
the current Sun quite well, to within 10% (see Fig. |8|. The 
remaining discrepancy is due to the slight variation in back- 
ground density, and the magnitude of the equilibrium tide, 
over the launching region. We note that the main change 
in the energy flux from modifying the profile of N'^ in the 
transition region is to change its frequency dependence. 

As can be seen from Fig. |8] the numerically computed 
asymptotic value of P differs from Fx by a factor < 2, even 

^ However, no simple physical arguments for such a profile have 
been found, which is perhaps to be expected since stellar models 
contain complicated combinations of physics. 
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in the case of short-period forcing. We would expect Fx to 
be correct to within a factor ^ 2 even if the slope varies by 
an order of magnitude, since it is only raised to the power 
— |. However, our slight refinement improves this estimate 
by - 50% for P = 1 d. 

We have performed numerical integrations for several 
different forcing frequencies (planetary orbital periods), for 
which we find F oc p~^-^^ for fixed stellar and planetary 
properties (other than the orbital period). If we take into 
account the fact that dx^r ^ ^ ^ P~^i and then consider 
a fixed tidal potential, we find that F oc cj^ *^^, which is 
slightly different from the power law dependence in Eq. [5l] 
The discrepancy most likely results from the variation in the 
background density, and the magnitude of the equilibrium 
tide (which forces the dynamical tide), within this region. 

For P < 2 d, which is the relevant regime for which this 
process is potentially important for the survival of close-in 
planets (BOlO), there is a few per-cent variation in back- 
ground parameters over a lengthscale Liaunch- This means 
that our calculation of F^ differs from the numerically com- 
puted value of P, by an amount that increases as P is made 
smaller to a maximum of 20% when P = 0.5 d. Neverthe- 
less, this discrepancy is small, and our semi-analytic esti- 
mate F^ is a good approximation to P for all cases that 
we have modelled. This is used to provide an estimate of 
and thus the tidal torque, in § [7| 



6.4 Variation between different stellar models 

Eq. [5l] matches the numerically computed value to within a 
few per-cent for a variety of solar-type stars. This is because 
it generally arises that N'^ oc x", with a ^ 0.5, when the 
launching region is a few percent of the stellar radius. This 
occurs for the waves excited by planets in short-period or- 
bits. We have confirmed this by computing P in a number of 
stellar models with masses in the range 0.5 ^ ui^/Mq ^1.1, 
and ages that represent the range of main-sequence ages ex- 
pected for these stars. These were computed using ASTEC 
( Christensen-Dalsgaard||2008| . 

The main uncertainty in these models is the profile of 



N within the launching region, especially since the slope 
within O.O2P0 of the interface is not well constrained by 
theory or observations. Helioseismic observations are not yet 
able to constrain the stratification within this region, owing 
to the lack of observed g-modes { Ellis 1984"; 'Appourchaux 
et al. 2010). In addition, the relevant physics included in 



the stellar models is also uncertain, particularly as a result 
of changes to the compositional gradient from convective 
overshoot. The inclusion of helium settling tends to make 
the interface profile sharper, and the inclusion of turbulent 
diffusion, that results from convective overshoot, and often 
parameterised using a simplified ID model of this process, 
tends to smooth out the profile near the interface (J0rgen 
Christensen-Dalsgaard, private communication). However, 
these changes are small and occur only within a region 
smaller than the launching region for P < 3 d. This means 
that uncertainties in the observations, and the physics, at 
the interface are unlikely to significantly change P for P < 3 
d, which are the planets whose survival could be threatened 
by our mechanism (as we discuss below). In addition, |^^| 
is only raised to the -2/5 power in our model. 
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Figure 10. In the top panel we plot N"^ in a I.OM0 star at 
t = 6.95 Gyr, with our fitted solutions in the bottom panel. This 
model has a "bump" in the N"^ profile, but this occurs over a 
region smaller than Liaunch^ cind so does not reduce the accuracy 
of our analytic model (from Eq. |B16| > in the launching region (or 
the corresponding energy flux), shown in the bottom panel. 



For Model S, the compositional gradient (V^) is unim- 
portant in the launching region, and N'^ primarily results 
from temperature gradients. As the star evolves, the set- 
tling of elements heavier than H, produces a compositional 
gradient at the top of the RZ, which can produce "bumps" 
in the profile of N'^ . In Fig. 



10 



we plot an example of the 
profile of N'^ for a IM© model (with Z = 0.02) at t = 6.95 
Gyr, together with the best fit solution in the launching re- 
gion. For P < 3 d, Liaunch is generally much larger than 
these "bumps" , so the wave launching process does not no- 
tice such departures from a smooth stratification profile, and 
Eq. |5l] remains a good approximation for the energy flux. 
However, if the frequency is sufficiently low (P > 6 d), the 
radial extent of the launching region can become comparable 
with the size of these "bumps" , and the numerical solution 
can depart appreciably from our analytical model. Planets 
in such orbits are very unlikely to be affected by stellar tides 
at such orbital distances, so we do not consider such effects 
worthy of further consideration. To summarise, we have con- 
firmed that Eq. [5l] is a good approximation for the energy 
fiux for planets on orbits of a few days throughout the range 
of solar-type stars in our study. 



6.5 Are ID linear hydrodynamic calculations 
reasonable? 

The calculations done in § |6] are performed under the as- 
sumption of linearity. We have demonstrated in § |5] and in 
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BO 10, that this is not a vahd assumption near the centre of 
a star. A simple estimate of the nonhnearity in IGWs in the 
launching region, for want of a better measure, compares 
the radial displacement to the radial wavelength, for which 
^t/Xr ^ 10~^, for a hot Jupiter orbiting the Sun on a one- 
day orbit (this estimate can be obtained from Fig.|7|. This 
value increases linearly with the mass of the planet, but we 
are still in the linear regime even if we have a close-binary 
perturber, which indicates that linearity is likely to be a 
good approximation in the launching region (and through- 
out the RZ, except for the central regions). 

However, it remains to be seen whether F for the ID 
calculations is the same as that in 2D or 3D numerical sim- 
ulations of realistic tidal forcing in a model including both 
a CZ and a RZ. Such simulations as Rogers et al. (2006) 



could be performed of the whole star subject to tidal forcing. 
The turbulent convection in these simulations will produce a 
spectrum of waves at the top of the CZ, in addition to those 
excited by tidal forcing. It may be that the interaction of 
these waves reduces the ingoing energy flux. Alternatively, 
the profile of in the transition region could be modi- 
fied by realistic modelling of convective overshoot, altering 
the strength of the stratification within a few percent of a 
pressure scale height from the interface. This would affect 
F if the overshoot region is comparable with the size of the 
launching region, though probably not significantly, as we 
discussed in the previous section. 

A toroidal magnetic field in the launching region could 
also affect the amplitudes of these waves, and the value of 
F. [Rogers MacG regorj ([Mq]) find that when the IGW 
frequency is approximately equal to the Alfven frequency 
(cja), strong wave refiection occurs. This behaviour follows 
from the dispersion relation for IGWs in the presence of 

11999'), and could have 



magnetic field (e.g. Kumar et al 



important consequences for F, if the magnetic field is suffi- 
ciently strong. However, for uj ^ uja the launching region, 
we require = 27r y^JIopbrb / P > 6 MG, when P = 1 d, 
which is close to being ruled out on empirical grounds, from 
measurements of the solar oblateness (Triedland & Gruzinov 
[2004 ). U (jJa < (jO, little attenuation of wave energy over the 



non-magnetic case is found by Rogers & MacGregor (2010), 



so this seems unlikely to affect F in our case. 

In the innermost wavelength, a strong magnetic field 
would be able to refiect IGWs before they reach the centre, 
if cj ~ CJA in this region. This requires a toroidal field of 
strength > 6 MG, or a poloidal (radial) field of strength 
Br > 28 MG. Since it is unlikely that such fields could exist 



in the RZ (Friedland & Gruzinov 2004), a magnetic field 



will probably not affect the refiection of IGWs (excited by 
short-period planets) in the RZ. It is therefore appropriate 
to ask whether the waves will break on reaching the centre. 



7 TIDAL Q' AND THE BREAKING 
CRITERION 

The modified tidal quality factor of the star that results 
from critical layer absorption can be obtained from Eq. I 



51 



(a 

slight modification from BO 10, resulting from the adopted 



profile of N"^ in the launching region), evaluating to 



^9x 10^ 



1 day 



(52) 



with the exact value depending on the stellar model adopted. 
However, this value is found to vary by only a factor of 5 
throughout the range of main-sequence stars in our study, 
for a given orbit. As the star evolves on the main-sequence, 
the position of the interface rt moves inwards towards higher 
density material, which slightly increases F. This means that 
Q' tends to decrease with main-sequence age, though only 
by up to a factor of 5. The dissipation that results from wave 
absorption at a critical layer thus becomes more effective as 
the star evolves. 

The waves break if 







C 



Mr. 



> 3.6, 



(53) 



Cq J \Mj J J \1 day ^ 

where 1/2 ^ /3 ^ 2 is a coefficient]^ which is unity for 
the current Sun but varies with the stellar model, and 
Cq = 8 X 10~^^m~^s~^. This expression is valid if we are 
sufficiently far from resonance with a global mode. Note, 
however, that if we are close to a resonance, the central 
amplitude may be much larger than these estimates would 
predict, which would make breaking more likely. Note also 
that this criterion becomes much easier to satisfy as the 
star evolves, since C increases with main-sequence age (see 
Fig. 1 in BOlO). If a planet exceeds this criterion, then its 
orbital migration is rapid (^ Myr; see Eq. 105 in BO 10) 
once critical layer absorption occurs, and would therefore 
threaten the survival of any short-period planet that causes 
wave breaking with P < 3 d. Very massive perturbers, with 
larger orbital moments of inertia than the spin moment of 
inertia of the RZ of the star, may be able to synchonise the 
spin of the RZ with the orbit, and prevent further inward 
migration. This may be an explanation for the synchronisa- 
tion of close-binary stars (e.g. Mazeh|2008 ). 



7.1 The survival of short-period planets 

Our most important result is that this mechanism can po- 
tentially explain the survival of all short-period extrasolar 
planets around solar- type star^ with masses in the range 
0.5 ^ ui^/Mq ^1.1. From using the closest fit stellar model 
to each of these stars, we find that no planet clearly sat- 



isfies Eq. 53 with a period P < 3 d (planets much further 
out may satisfy the criterion, though this tidal effect is then 
unimportant). All of these planets are insufficently massive 
(or orbit sufficiently young stars, with low values of C /Cq) 
that they are unlikely to cause wave breaking at the centre 
of their hosts. The dominant mechanism of tidal dissipa- 
tion in these stars is therefore likely to be damping of the 
equilibrium tide by turbulent convection. This is likely to 
be relatively inefficient for these periods, since the turbulent 
viscosity must be reduced when the orbital period is shorter 
than the convective timescale (see e.g. Zahn|1966 Goldreich 



^ Which was written as (g/gp)^/^ in BOlO Eq. 103; see that 
paper for further details. 



' Contained in the catalogue at |jittp://www.exoplanet.hanno- 
rein.de/ or |http://exoplanet.eu/catalog.php 
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k Nicholson|1977| [Goodman Oh|1997||Penev et al.|2007 ). 

This could be an important explanation for the survival of 
these planets. 

This mechanism can also explain why the most massive 
short-period planets (which still have smaller or compara- 
ble moments of inertia in the orbit as the RZ of the star), 
such as WASP- 18 b ( [Hellier et alT2 Q09), WASP- 14 b ( |joshi 
|et a l. 2009) or CoRoT-14 b, are exclusively found around 
F-stars, which have convective cores, and in which critical 
layer formation induced by wave breaking at the centre is 
unable to operate. In addition, these stars have been found 
to have weaker tidal dissipation in their outer CZs than G- 
stars (Barke r Ogilvie||2009p . Note, however, that planets 
with much larger masses may have sufficient orbital angular 
momentum to be able to synchronise their stars and reach an 
approximate tidal equilibrium state, neglecting stellar mag- 



Schubert|1967||Knobloch Spruit[T982t , or magnetic in- 
stabilities, such as the magnetorotational instability ( [Balbus 
&: Hawley|1994). However, these mechanisms are unlikely to 



netic braking, which would prevent orbital decay (e.g. Lev- 
rard et ar] |2009), even for planets in orbit around a solar- 
type star. Consequently, we make the prediction that fewer 
massive planets in the range 3 < Mj < 20, in orbits with 
P < 2 — 3 d, that satisfy Eq.|53] will be found around solar- 
type stars. 



8 CONCLUSIONS 

In this paper, we have presented the results of a set of 3D 
simulations of IGWs approaching the centre of a solar-type 
star. The aim of this work is to determine the importance 
of a nonlinear mechanism of tidal dissipation in solar-type 
stars, first proposed by GD98, continuing the investigation of 
BO 10. We have confirmed that the main results of BO 10 are 
unaffected by the extension to 3D by performing numerical 
hydrodynamical simulations, using a Boussinesq-type model 
solved with a pseudospectral method. 

We first derived a linear wave solution in 3D, and found 
that nonlinearities do not vanish for this wave, unlike the 2D 
solution, which is exact. Nevertheless, these waves are found 
to reflect approximately perfectly for moderate-amplitudes, 
a result which we have qualitatively confirmed in numerical 
simulations of moderate-amplitude tidal forcing. 

The general picture for high-amplitude forcing is that 
IGWs break within the innermost wavelengths of a star, if 
they reach the centre with sufficient amplitude to overturn 
the stratification. If this occurs, they form a critical layer, 
which we have confirmed from the simulations efficiently ab- 
sorbs ingoing wave angular momentum. This results in the 
star being spun up to the orbital angular frequency of the 
planet, from the inside out. This could be very important to 
the survival of massive planets in short-period orbits around 
solar- type stars. 

One noticeable difference in 3D is that the absorption of 
/ = m = 2 ingoing waves results in the formation of latitu- 
dinal differential rotation. This is perpetually reinforced by 
critical layer absorption. Instabilities may act on this rota- 
tion profile, which could homogenise the horizontal angular 
momentum distribution. These include shear instabilities. 



which can be linear (Watson 1981) or nonlinear instabil 



be able to prevent the critical layer absorption, and thus 
prevent the tidal engulfment of a short-period planet. 

In these simulations we artificially forced waves at the 
top of the computational domain. In reality, the waves are 
launched at the top of the RZ. In this paper we also solved 
the linearised equations for the adiabatic tidal response, in 
order to study the infiuence of the launching region in more 
detail. We provide an expression for the energy fiux Eq. |5H 
which has been checked to agree very well with the numeri- 
cally computed energy fiux for all orbits for which this effect 
is potentially important, and in a wide range of solar- type 
stellar models. We then presented the breaking criterion that 
must be satisfied for the waves to break, slightly refining 
BOIO. 

Finally, in §[7|we presented the tidal Q' that results from 
this process, together with a breaking criterion that deter- 
mines when these waves break. This allowed us to outline 
a possible explanation for the survival of all short-period 
extrasolar planets, that is not in confiict with current ob- 
servations, and makes predictions which will be tested by 
ongoing observational studies such as WASP and Kepler. 

As in BOIO, our simulations do not show any instabil- 
ities that act on the waves when they have insufficient am- 
plitude to overturn the stratification. This result suggests 
that to obtain any efficient tidal dissipation, the waves must 
have sufficient amplitude to satisfy Eq.[53] so that they over- 
turn the stratification. However, it may be that weaker para- 
metric instabilities operate for waves with lower amplitudes 
(suggested by GD98). A detailed stability analysis of the 2D 
exact wave solution written down in BOIO is ongoing, and 
should shed light on this matter. 

In our explanation in the previous section, we ignored 
the potential effects of the evolution of either the stellar 
eigenfrequencies or the tidal frequency, that could result in 
a passage through a resonance with a global mode of oscil- 
lation. If this occurs, and simple estimates suggest that this 
is very likely at least once in the lifetime of a system, then 
the waves may break, even if the amplitude of the tide is 
small before the passage through resonance. Future work is 
required to study this problem in detail. 

The simulations performed in this paper used a 
Boussinesq-type model, which is only valid in the central 
< 5% by radius of a solar-type star. For future work, it 
is possible to perform simulations using an anelastic code, 
which could take into account the radial variation in stellar 
properties throughout the whole RZ (e.g. [Rogers &: Glatz^ 
maier|2005 ) . We have also neglected the rotation (or possible 



ities, that set in at a critical Reynolds number (Richard 



&: Zahn[l999 ). These have growth times comparable to the 
tidal period, and could transfer angular momentum latitudi- 



nally. There are also doubly diffusive instabilities ( Goldreich 



differential rotation) of the RZ, which is reasonable since we 
are considering the waves excited by short-period planets 
in slowly rotating stars. However, the inclusion of rotation 
(and differential rotation) would certainly be interesting to 
study for more rapidly rotating stars. One consequence could 
be that the latitudinal differential rotation observed in the 
simulations in 3D may be enhanced by rotation, because 
inertia-gravity waves can be confined in the equatorial re- 
gions. 

In addition, we have so far considered a planet on a cir- 
cular, coplanar orbit. Since many planets have been observed 
with eccentric or inclined orbits, this makes a study of the 
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circularisation and potential spin-orbit alignment through 
wave breaking and the resulting critical layer absorption, 
seem an attractive avenue of future research. 
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APPENDIX A: IW/OW DECOMPOSITION 

In § |3] we derived an analytic solution for the waves near 
the centre. We can use this to deconstruct the numerical 
solution into a single IW and OW, as in BO 10. Doing this 
enables us to quantify the amount of angular momentum 
absorbed as these waves approach and/or reflect from the 
centre. However, the reduction in resolution per dimension 
required for the transition to three dimensions means that 
an alternative method is preferable, which includes the in- 
formation at several points, to reduce numerical error. The 
approach we use is now described. 

At every point with spatial coordinates (x, y, z), we have 
four pieces of information, namely, Ur, ue, U(f) and b. Thus, 
it is possible to compute the (complex) IW and OW am- 
plitudes Ain and Aout for the I — m — 2 wave at a single 
point. This is, however, computationally expensive, since the 
routines for computing Bessel functions are relatively slow 
(typically several hundred times slower than a square root). 
In addition, this method would be subject to potentially 
significant round-off errors. 

We fit the simulation output to a linear model, cor- 
responding to a an IW and an OW. This problem is an 
overdetermined system of linear equations in terms of the 
unknown wave amplitudes, which we write as 



(Al) 



where y — {ur,ue,Uci),b) is a vector of data variables at 
each grid point, of size 4A^, where N — NxNyNz is the 
total number of grid points. M is the number of spherical 
harmonics for which we compute the wave amplitudes, which 
is usually taken to be one. x = {Ain, Aout) is a vector of size 
4M < 4A/', whose components are the (complex) IW/OW 
amplitudes (for each / and m value) , whose values are to be 
determined. The matrix A contains the IW and OW radial 
functions and spherical harmonic functions, evaluated at the 
selected grid points using the linear solution in § [3] and has 
size 4iV X 4M. 

We use a method of least squares to fit our model to 
the data (Press et al.| |1992| . This finds the best fit between 
the linear model data and simulation data, and computes the 
solution for which the sum of the squared residuals is its least 
value. This problem has a unique solution if the AN columns 
of the matrix A are linearly independent, which is true in our 
case, since 

Vr. The solution is found by solving the normal equations 



We take into account radial variations in the amplitudes by 
splitting up the region inside the forcing region into a set 
of concentric spherical shells of thickness Sr ^ Ar/2, after 
removing an inner region of a few grid cells. This approach 
assumes the solution is locally independent of r, hence we 
can ignore radial derivatives of the amplitudes within each 
shell. This is not valid in regions where the solution varies 
rapidly. In addition, we speed up computation by stepping 
over the grid points in each Cartesian direction, by factors 
istep, jstep, kstep, choscu to take values between 1 and 10. We 
always ensure that sufficient grid points are available in each 
shell to accurately compute the amplitudes. 

A Matlab routine which reads in SNOOPY data and 
solves the linear least squares problem has been written, and 
tested with various combinations of analytic IW/OW solu- 
tions. We compute the reflection coefficient IZ = Aout /Ain- 
For perfect standing waves, Ain — Aout, and IZ — 1. If the 
IW is entirely absorbed at the centre, then 7^ = 0. 

The main disadvantage of our approach is that uj ^ 
1,/ 7^ 2,m 7^ 2 components also contribute to the ampli- 
tudes. This problem can be ignored if we trust the computed 
values of IZ only where the solution is well described by the 
linear solution, i.e., far from the wave breaking and forcing 
regions. We therefore reconstruct Ur, and plot it along the 
X-axis, ue along the line y — U(j) along the line y — 
and h along the line y — x. All components of the solution 
are not zero for all radii along these lines when ^ = (this 
requires shifting the phase of the linear model solution when 
the wave phase is nonzero in the simulation data) . The sim- 
ulation variables and IW/OW solutions are plotted using 
this method in Figs. ^ and ^ Comparing the reconstructed 
solution with the simulation data allows us to check whether 
the deconstruct ion has worked. 



APPENDIX B: ANALYTIC CALCULATION OF 
F IN THE LAUNCHING REGION 

In this section, we adopt the notation of OL04, to avoid 
reproducing the results in that paper, since we are simply 
extending the results of their §4.4. In the launching region, 
we want to match the solutions of the RZ and CZ. Near 
the boundary r = r^, we assume N'^ oc (r^ — r)", where 
a G M^, for the moment left unspecified. The characteristic 
radial extent of the transition region is of order e^, where 
P G M^, and e is the ratio of the spin frequency to the 
dynamical frequency, which is much smaller than unity for 
slowly rotating bodies. We write 



r — rh — ex. 



(Bl) 



where x is an inner variable, which is of order unity in the 
launching region. We can write 



(A2) 



2olI3\ 



(B2) 
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where V = 4^ > 0. If we pose the perturbation expansions 



/ 




(B3) 


u'e 


~ eu0{x,6), 


(B4) 






(B5) 


1 

P 




(B6) 


1 

V 




(B7) 






(B8) 



then we obtain the hnearised system of equations 111-116 
from OL04, except that we replace x by x^ in the buoyancy 
equation (Eq. 115). A Hnearised equation for the modified 



pressure {W = ^ 
now be derived, for which 

+ (x-'^d^w) = 0, 



+ in the launching region can 



(B9) 



where all coefficients are evaluated at r = r?,. This equation 
can be solved using the separation of variables 



(BIO) 



The operator C contains only angular derivatives, and 



Cwi — XiWi, where 



and k is the WKB wavenum- 



ber of the waves. We thus obtain the equation 

+ =0, (Bll) 

where we have defined z — kix, and the length scale 



VXi 



(B12) 



In the text this lengthscale is referred to as Liaunch- The 
solution of this equation can be written down as a linear 
combination of Bessel functions of the first and second kinds, 
of order 

Oi-\-2 



aiZ "^2^ \ J g+i ( ^ ^ 2:^2" I + 



^+2 \2 -\- a 



1 L (B13) 



^ \2 + a 

so that the complete solution in the transition region is 



W 



, 2 2+0 
Sji Ya+i I — z 



} 



a+2 \Z -\- a ' ' 

This solution should match onto the WKB solution in the 
RZ when x ^ 1, constraining /3 as a function of a. This 
arises from considering the asymptotic form of the phase 
at large whose radial derivative should match the WKB 
wavenumber in the RZ. From this we find /S = 2/(2 + a). 

If we from now on restrict ourselves to a = 1/2, then 
the corresponding complete solution in the transition region 
is the sum 

W = J2aizi jjs +SiiJ3 (^^z^^^WiiO). (B15) 

Similarly we can write the radial displacement as 

= ^6,7^1^1 (^^z^^±itiJ2 (^^zi^'^WiiO). (B16) 



for a superposition of wave solutions in this model. 

Deep in the RZ, our solution should reduce to a wave 
with an inwardly directed group velocity, for which Si = ±1 
(sign depending on that of the frequency). The asymptotic 
forms of the radial part of this solution at large z is 



J, ( -.-^ 



(-1)- 



4 5 



'yj.iexp|^.t|.(B17) 



Matching of at r = requires the solution in the CZ to 
be continuous with the solution in the transition region at 
z = 0. This determinines the amplitudes a^. The inward en- 
ergy flux in (inertia-) gravity waves is found from evaluating 



r sin OdO 



(BI8) 



[r(i)] 



rPoCJ(i)Sgn(a)) 



r \ 
V 



i Jo \wi\^smede 

which can be compared with OL04 Eq. 124. We find in gen- 

8 + 3cx 

eral that F oc oj 2+« . Thus we see that the main change in 
the energy flux from modifying the profile of N'^ in the tran- 
sition region is to change its frequency dependence. There 
are also 0(1) changes to the numerical factors. 
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